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The relative dispersion of pairs of inertial particles in incompressible, homogeneous, and 
isotropic turbulence is studied by means of direct numerical simulations at two values of 
the Taylor-scale Reynolds number Re\ ^ 200 and Rex ~ 400, corresponding to resolu- 
tions of 512"^ and 2048'^ grid points, respectively. The evolution of both heavy and light 
particle pairs is analysed at varying the particle Stokes number and the fluid-to-particle 
density ratio. For particles much heavier than the fluid, the range of available Stokes 
numbers is St S [0.1 : 70], while for light particles the Stokes numbers span the range 
St £ [0.1 :3] and the density ratio is varied up to the limit of vanishing particle density. 
For heavy particles, it is found that turbulent dispersion is schematically governed by 
two temporal regimes. The first is dominated by the presence, at large Stokes numbers, 
of small-scale caustics in the particle velocity statistics, and it lasts until heavy particle 
velocities have relaxed towards the underlying flow velocities. At such large scales, a 
second regime starts where heavy particles separate as tracers particles would do. As a 
consequence, at increasing inertia, a larger transient stage is observed, and the Richard- 
son diffusion of simple tracers is recovered only at large times and large scales. These 
features also arise from a statistical closure of the equation of motion for heavy particle 
separation that is proposed, and which is supported by the numerical results. 
In the case of light particles with high density ratios, strong small-scale clustering leads 
to a considerable fraction of pairs that do not separate at all, although the mean separa- 
tion increases with time. This effect strongly alters the shape of the probability density 
function of light particle separations. 



1. Introduction 

Suspensions of dust, droplets, bubbles, and other finite-size particles advected by in- 
compressible turbulent flows are commonly encountered in many natural phenomena (see, 
e.g., |Csanady 1980|lEat"on fc Fessler 1994[|Falkovich et al. 20021IPost fc Abraham 2002[ 



IShaw 20031 IToschi fc Bodenschatz 2009p . Understanding their statistical properties is 
thus of primary importance. From a theoretical point of view, the problem is more com- 
plicated than in the case of fluid tracers, i.e. point-like particles with the same density as 
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the carrier fluid. Indeed, when the suspended particles have a finite size and a density ra- 
tio difi^erent from that of the fluid, they have inertia and do not foUow exactly the flow. As 
a consequence, correlations between particle positions and structures of the underlying 
flow appear. It is for instance well known that heavy particles are expelled from vortical 
structures, while light particles tend to concentrate in their cores. This results in the 
formation of strong inhomogeneities in the particle spatial distribution, an effect often 
refered to as preferential coneentration (see |Douady et al. 1991[ [Squires fc Eaton 199H 
lEaton fc Fessler 1994[) . This phenomenon has gathered much attention, as it is revealed 
by the amount of recently published theoretical w ork ([Balkovsky et al. 2001 IZaichik et al. 20031 
IFalkovich fc Pum ir 2004^ , and numerical studies (IColhns fc Keswani 20041 IChun et al. 20051 
iBec et al. 2007, Goto fc Vassilicos 2008[) . Progresses in the statistical characterization of 
particle aggregates have been achieved by studying particles evolving in stochastic flows 
by lSigurgeirsson fc Stuart 2002HMeh hg fc Wilkinson 2004|[B ec et al. 2005 l lOUa 2002l and 
in two-dimensional turbulent flows bv Boffctta et al. 20041 Also, single trajectory statis- 
tics have been addressed both numerically and experimentally for small heavy particles 
(see, e.g.. lBec et al. 200 6','Cenc ini et al. 20061 [Gylfason et al. 20061 IGera shchcn ko et al. 20081 
[Zaichik fc Alipchenkov 2008,,Ayyalasomayajula et al. 2008 IVolk et al. 2008,) . and for large 
particles ( Qureshi et al. 2007 IXu fc Bodenschatz 2008P . The reader is refered to lToschi fc Bodenschatz 20091 
for a review. 

In this paper we are concerned with particle pair dispersion, that is with the statistics, 
as a function of time, of the separation distance R{t) = Xi{t)—X2{t) between two inertial 
particles, labelled by the subscripts 1 and 2 (see IBec et al. 20081 IFouxon fc Horvai 20081 
IDerevich 20081 for recent studies on that problem). In homogeneous turbulence, it is 
sufl&cient to consider the statistics of the instantaneous separation of the positions of the 
two particles. These are organised in different families according to the values of their 
Stokes number St, and of their density mismatch with the fluid, f3. 

For our purposes, the motion of particle pairs, with given {St, (3) values and with 
initial separations inside a given spherical shell, R = \Xi{to) — X2{to)\ € [Rq, Rq + dRg] 
is followed until particle separation reaches the large scale of the flow. With respect to 
the case of simple tracers, the time evolution of the inertial particle pair separation R{t) 
becomes a function not only of the initial distance Rq, and of the Reynolds number of 
the flow, but also of the inertia parameters (St, (3). 

A key question that naturally arises is how to choose the initial spatial and velocity 
distributions of inertial pairs. Indeed, it is known that heavy (resp. light) particles tend 
to concentrate preferentially in hyperbolic (resp. elliptic) regions of the advecting flow, 
with spatial correlation effects that may extend up to the inertial range of scales, as 
shown in lBec et al. 200T1 Moreover, when inertia is high enough, the particle pair velocity 
difference, SuV = \Vi{Xi{t),t) — 'V2{X2{t),t)\, may not go smoothly to zero when the 
particle separations decreases, a phenomenon connected to the formation of caustics, 
see [Wilkinson fc Mehlig 20051 IFalkovich fc Pumir 20071 In our numerical simulations, 
particles of different inertia are injected into the flow and let evolve until they reach a 
stationary statistics for both spatial and velocity distributions. Only after this transient 
time, pairs of particles with fixed intial separation are selected and then followed in the 
spatial domain to study relative dispersion. 

By reason of the previous considerations, the main issue is to understand the role 
played by the spatial inhomogeneities of the inertial particle concentration field and by 
the presence of caustics on the pair separations, at changing the degree of inertia. We 
remark that these two effects can be treated as independent only in the limit of very 
small and very large inertia. In the former case, particles tend to behave like tracers and 
move with the underlying fluid velocity: preferential concentration may affect only their 
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separation. In the opposite limit, particles distribute almost homogeneously in the flow: 
however, due to their ballistic motion, they can reach nearby positions with very different 
velocities (IFalkovich et ai. 2002|) . In any other case of intermediate inertia, both these 
effects are present and may play a role in the statistics of inertial pair separation. 
It is worth anticipating the two main results of this study: 

(i) The separation between heavy particles can be described in terms of two time regimes: 
a first regime is dominated by inertia effects, and considerable deviations from the tracers 
case arise in the inertial relative dispersion; in the second one, the tracers behaviour is 
recovered since inertia is weak and appears only in subdominant corrections that vanish 
as 1/t. The crossover between these two regimes defines a new characteristic spatial and 
temporal scale, connected to both the size of caustics and the Stokes number, which 
influences the particle separation for not too long time-lags and not too large scale. 

(ii) The strong clustering properties that are typical of light particles may lead to the fact 
that many pairs do not separate at all: their statistical weight is clear in the separation 
probability density function (PDF), which develops a well defined power-law left tail. 

It would clearly be also interesting to investigate the dependence upon the Reynolds 
number of the inertial particle pair separation. Small-scale clustering seems to be poorly 
dependent on the degree of turbulence of the carrier flow (jCoUins fc Keswani 2004[ 
IBec et al. 2007|) . while much less is known about the Reynolds number dependence of 
the caustics statistics. Our numerical data do not allow to explore this question in detail, 
so that we will restrict ourselves to show data associated to the two Reynolds numbers 
in all cases when differences are not significative. 

In the case of fluid tracers, the standard observables are the time evolutions of the 
mean square separation and of the separation probability density function, for which 
well established predictions exist since the pioneering work of lRichardson 19291 We con- 
trast these observables obtained for tracers with the results for heavy and light inertial 
particles. 

The paper is organised as follows. In Sj2l we briefly recall the basic equations of motion 
and describe the numerical simulations. In ij3l we analyse the stationary distribution of 
heavy particle velocity differences, conditioned on the particle initial separation, high- 
lighting both the presence of small-scale caustics and the effects of particle inertia at 
those scales corresponding to the inertial range of turbulence. In ij4] we study the be- 
haviour of the mean separation distance of heavy pairs, at changing the Stokes number 
St; we also analyse the influence of the caustics in the initial statistics on the subse- 
quent pair separation evolution. A mean-field model, which is able to capture the main 
numerical findings, is proposed in the same section. The time evolution of the separation 
probability density functions is discussed in |j5]and we present the data for light particles 
in ij6l In ij7]we summarise the main findings. 

2. Equation of motion and numerical details 

We present results from direct numerical simulations of turbulent flows seeded with 
inertial particles. The flow phase is described by the Navier-Stokes equations for the 
velocity field u{x,t) 

dtU + u-Vu= -Vp + vV'^u + f, V-M = 0. (2.1) 

The statistically homogeneous and isotropic external forcing / injects energy in the first 
low wave number shells, by keeping constant their spectral content fsee lChen et al. 1993p . 
The kinematic viscosity v is chosen such that the Kolmogorov length scale rj 6x, where 
Sx is the grid spacing: this choice ensures a good resolution of the small-scale velocity 
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N 


Rex 


V 


5x e 


V 




^dump 


&t 


Run I 
Run II 


512 
2048 


185 
400 


0.01 
0.0026 


0.012 0.9 
0.003 0.88 


0.002 
0.00035 


0.047 
0.02 


0.004 
0.00115 


0.0004 
0.000115 



Table 1. Eulerian parametres for the two runs analysed here: Run I and Run II in the text. N 
is the number of grid points in each spatial direction; Rex is the Taylor-scale Reynolds number; 
r\ is the Kolmogorov dissipative scale; &x = C-jN vs, the grid spacing, with C = 2n denoting the 
physical size of the numerical domain; = \AV^ ^'he Kolmogorov dissipative time scale; e is 
the kinetic energy dissipation; v is the kinematic viscosity; Tdump is the time interval between 
two successive dumps along particle trajectories; 5t is the time step. 

dynamics. The numerical domain is cubic and 27r-periodic in the three directions of space. 
We use a fully dealiased pseudospectral algorithm with 2"^^ order Adam-Bashforth time- 
stepping (for details see lBec et al. 20061 ICencini et al. 2006p . We performed two series of 
DNS: Run I with numerical resolution of 512'^ grid points, and the Reynolds number at 
the Taylor scale Re\ « 200; Run II with 2048^ resolution and Re\ « 400. Details of the 
runs can be found in Table [1] 

The particle phase is constituted by millions of heavy and light particles — the latter 
only for Run I — with different intrinsic characteristics. Particles are assumed to be with 
size much smaller than the Kolmogorov scale of the flow, 77, and with a negligible Reynolds 
number relative to the particle size. In this limit, the equations ruling their dynamics 
take the particularly simple form: 

X = V, V = -—[V -u{X,t)]+ [iDtu{X,t) , (2.2) 

Ts 

where the dots denote time derivatives. The particle position and velocity are (-X'(t), V{t)), 
respectively; u{X{t),t) is the Eulerian fluid velocity evaluated at the particle position, 
and DfU is the so-called added mass term, which measures the fluid acceleration along 
particle trajectory. The adimensional constant [3 = 3pf/ {pf + 2pp) accounts for the added 
mass effect through the density contrast between particles pp and fluid pf. The parti- 
cle response time, appearing in the Stokes drag, is Tg = 2ppa^ / (9pfv), where a is the 
particle radius. Particle inertia is quantified by the Stokes number that is defined as 
St = Ts/Trj, where = (i^/e)^/^ is the flow Kolmogorov timescale and e the average 
rate of energy injection. Equation (|2.2p has been derived in |Maxey fc Riley 1983| under 
the assumption of very dilute suspensions, where particle-particle interactions (collisions) 
and hydrodynamic coupling to the flow can be neglected. 

For Run I, we show results for the following set of {St,P) families: (i) very heavy 
particles [/3 = 0]: St = 0.0,0.6,1.0,3.3; {ii) light particles [/3 = 2,3]: St = 0.3,1.2,4.1. 
For each family the typical number of particle pairs that are followed is around 5 x 10'''. 
For Run II, we show results only for heavy particles but with a larger range of variation in 
the Stokes number: St = 0.0, 0.6, 1.0, 3.0, 10, 30, 70. Typical number of particle pairs for 
each family is ~ 10^. Once injected particles have relaxed to their steady-state statistics, 
pairs have been selected with the following initial separations: Rq < rj and Rq S [4:6]r7 
for both Run I and Run II, and Rq G [9 : 11] 77 for Run II only. 

Beside the time evolution of particle pairs, we also have instantaneous snapshots of 
the two phases (fluid and dispersed), with a much higher particle statistics: around 10^ 
per family for Run I, and 10^ per family for Run II. These are used to measure the 
stationary — i.e. not along the trajectories — distribution of particle velocity increments 
discussed in next section. 
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3. Stationary distributions: velocity increments conditioned on 
particles separation 

Turbulent pair dispersion for tracers is classically based on the application of similarity 
theory for Eulerian velocity statistics: depending on the value of space and time scales, 
velocity increment statistics differently affect the way tracers separate. This results in 
different regimes for relative dispersion, see e.g. ISawford 20011 

In the case of inertial particles, the same reasoning holds, so that to analyse the way 
inertial pairs separate in time, the stationary statistics of particles velocity differences 
has to be investigated first. A stationary distribution for the typical velocity differences 
between two inertial particles is obtained by imposing periodic boundary conditions inside 
the physical volume and then measuring velocities on such a thermalised configuration. 

We are interested in the scaling behaviour of velocity increments at varying the degree 
of inertia and the distance between the particles (in the dissipative or inertial range of 
the turbulent fluid flow). To fix the notation, we denote by Uq the typical large-scale 
velocity of the fluid tracers and by L the integral scale of the flow. Moreover we define 

iRYst = \Vi{Xi{t))-V2{X2{t))l (3.1) 

as the velocity difference at scale i?, conditioned on the presence of a pair of particles with 
Stokes number St, separated with a distance R — |JCi(io) — -X^2(^o)|- Since we arc here 
interested in the case of heavy particles only, the Stokes number is sufficient to identify 
a given particle family. For convenience, we introduce a specific notation for the tracer 
stationary velocity statistics: 5jiu = (5_RV(5t=o)j which is exactly equal to the Eulerian 
velocity increment at scale R. 

Recently, IBec et al. 20081 have shown that to describe inertial pair dispersion in syn- 
thetic flows it is useful to introduce the local or scale- dependent Stokes number, using the 
ratio between the particle response time and the typical eddy turnover time = R/5ru 
of the underlying fluid at a given scale: St{R) ~ Ts/tpi ^ TgSf/u/R. For real turbulent 
flow where different scaling ranges arc present, we can equivalcntly define a scale de- 
pendent Stokes number St{R) that recovers the usual definition of the Stokes number 
St{R) ~ 5t = Ts/tjj when i? <C r/ and behaves as St{R) ~ Tse'^-Z^R"'^/^ when R:$> rj. The 
typical behaviour of St{R) is sketched in Fig. [1] for two different values of the Stokes 
number St = 3, 70 and using a Batchclor-like parametrisation of the fluid velocity (see 
IMeneveau 1996^ : 

R 

^R"" = ^0 (^2 + j^2) 1/3 - (3.2) 

For Stokes numbers, St, order unity or larger, there always exists a typical scale where 
the local Stokes number, St{R), becomes order unity, 

R*{St)^r^St^'^. (3.3) 

Such a scale, which is well in the inertial range if the Stokes number St is sufficiently large, 
can be considered a rough estimate of the upper bound for the region of scales where in- 
ertia plays an important role in the particle dynamics. We expect that two main features 
might be important in characterising the inertial particle stationary velocity statistics 
SrV, with respect to that of tracers Sru. The first concerns the small-scale behaviour of 
the particle velocity statistics. At small scales R <^ rj and for large-enough Stokes num- 
bers, the presence of caustics makes the particle velocity increments not differentiable. 
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Figure 1. Left panels. Bottom figure: befiaviour of tfie scale-dependent Stokes number, St{R) 
as a function of the scale R normalised with the Kolmogorov scale 77, for two Stokes numbers 
St = 3,70 (bottom and top, respectively). The horizontal thick line is for St{R) = 1. Top 
figure; the scaling behaviour for the fluid tracer velocity increments versus the scale as given 
by (|3.2|l . Notice that the scales -R* where St{R*) = 1 fall in the inertial range of the Eulerian 
fluid velocity. Right panel: the function 'yiSt) defining the small scales power law behaviour of 
caustics statistics at changing inertia. Notice that for small values of the Stokes number St, 
7 — > 1, i.e. particle velocity is differentiable; at high inertia, 7^0 indicating the existence of 
discontinuities in the particle velocity increment statistics. 

This feature can be accounted for by saying that 

/ n\ 7(-St) 

SRVst-V^A-] ; i?«r;, (3.4) 

where the V^^ is a constant prefactor and the function ^{St) gives the typical scaling 
of caustic-like velocity increments. Indeed we do expect that at changing the inertia of 
the particles, the statistical weight of caustics might monotonically vary as follows: at 
small St, \imst^oj{St) = 1 , i.e. the value for smooth, differentiable Eulerian statistics 
of tracers; at large values St — > 00, it should approach the discontinuous limit j{St) — > 0, 
valid for particles that do not feci underlying fluid fluctuations at all. The right panel 
of Fig. [1] shows the typical shape of the function jiSt) that is expected to be valid for 
turbulent flows. 

The second important feature concerns the particle velocity statistics at scales larger 
than the scale R*{St) previously defined, but smaller than the integral scale of the fluid 
flow. For any fixed Stokes number and for a large-enough Reynolds number, we expect 
that inertia becomes weaker and weaker, by going to larger and larger scales R ^ R*{St). 
In such a case, particle velocity increments are expected to approach the underlying fiuid 
velocity increments: 

/ n\ 1/3 

SRVst^V^,SRU^V§,Uoi-] ; i?*(5i) « i? « L, (3.5) 

where for simplicity we have neglected possible intermittent correction to the Kolmogorov 
1941 (K41) scaling of the fluid velocity (see lFrisch 1995l for details). Clearly, the Reynolds 
number has to be sufficiently large to provide a well-developed scaling region R*{St) <C 
R<^ L, before approaching the large scale L. We emphasise that in (|3.5p . an adimensional 
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R/Tl R/Tl 

Figure 2. Right figure: particle velocity structure function of order p = 1 versus the scale R/rj, 
for various Stokes numbers, St = 0, 0.6, 1, 10, 70, and for Reynolds number Rex ~ 400, Run II. 
The statistics for fluid tracers {St = 0) correspond to the solid line. Statistical errors are of the 
order of twice the size of symbols for scales smaller then 77 and become comparable with the size 
of symbols in the inertial range of fluid velocity statistics. The differentiable scaling behaviour 
DC 7? in the dissipative range, and the Kolmogorov 1941 behaviour oc R^^^ in the inertial range 
of scales are also shown. Left figure: zoom-in of the inertial range, same symbols as right. Inset: 
behaviour of the amplitude prefactor, Vst as a function of the Stokes number St, as measured 
from the velocity increments at the integral scale L, Run II. 



normalisation factor Vg^. has been introduced: it takes into account possible filtering ef- 
fects induced by inertia at large scales. The normalisation is such that V^5t^o) ~ ^' '^hil^ 
for any Stokes larger than zero V^gt) — ^■ 

In Fig. [2] we test the validity of the previous picture by analyzing the typical velocity 
fluctuation, (|(5_RV5t|), at changing Stokes number and for data of Run II at Reynolds 
number Rex ~ 400. At small scales one detect the presence of caustics in the velocity 
statistics, vifith a non smooth scaling behaviour below the Kolmogorov scale 77. At scales 
within the inertial range and when the Stokes number is sufhciently large, the effect of 
caustics affects also particle velocity statistics, up to a characteristic scale which becomes 
larger and larger by increasing particle inertia. Beyond this scale, particle velocity incre- 
ments tend to approach the scaling behaviour of the fluid tracers, but their amplitude 
is depleted of a factor l/V^^, which increases with the Stokes number, as shown in the 
inset of the right-hand panel of Fig. [2l A similar behaviour is expected for higher order 
fluctuations, if wc neglect the role of intcrmittcncy. 

It is interesting to consider the scaling behaviour of particle velocity in terms of the 
underlying velocity statistics, not only at very small or very large separations, but for 
any value of the scale R. This is not straigthforward, since we have to account not only 
of the fluid Eulerian statistics at the dissipative and inertial range of scales, but also the 
modifications due to the inertia. This is responsible, as we have seen, for the appeareance 
of a new relevant scale, and for filtering effects in the velocity amplitude. 

To fully characterise particle velocity increments, we notice that the Stokes scale, 
R*{St) defines a typical Stokes-velocity: this is the fluid velocity increment at the Stokes 
scale, 6u*{St) ~ Sb*u (see left panel of Fig. [1]). Previous reasonings can be summarised 
in the following interpolation formula for the heavy particle velocity increment: 

S^Vst = Vl, (<5hu)^(^*(^» [{Snuf + ci {5u*{St)f] . (3.6) 

The above expression is a Batchelor-like paramctrisation but in the velocity space, with 
a transient velocity given by the Stokes velocity, 5u*{St). 
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Figure 3. Scaling behaviour of the particle velocity structure function of order one, versus the 
normalised scale R/i]. Solid lines: fit of the data of Fig. (2] Run II, using the interpolation formula 
l|3.6[l . Here the large scale prefactors Vgt are those measured on Run II of the simulation, and 
shown in the inset of the right panel of Fig. [2] Inset: enlargement of the crossover range, where 
Sru ^ SrV. 

Once known the large scale normalization function Vg^, the caustic exponent j{x) 
(introduced in lBec et al. 2005]) and the reference fluid velocity increment Sru, then the 
formula has one free parameter only. It is the prefactor ci appearing in front of the Stokes 
velocity Su*{St), whose value depends again on the inertia of the particles. 

In Fig. [21 we show the result of the fit in terms of the expression (|3.6p . where the caus- 
tics scaling exponent has been chosen as 7(2;) = [1 — 2/7ratan(j:)]: this functional form 
provides a good fit to the numerical results. Details of the small-scale caustic statistics 
will be reported elsewhere. 

The qualitative trend is very well captured by the interpolation function proposed. 
Notice that in (|3.6p . the argument of j{St) is not the simple Stokes number at the 
Kolmogorov scale, but the scale-dependent one St{R): j{St) j{St{R)). This further 
ingredient is needed to take into account the fact that in presence of a rough underlying 
fluid velocity, as it happens in the inertial range of scales, no simple power law behaviour 
is expected for the scaling of particle velocity statistics. This was previously remarked in 
IBec et al. 2008( in the study of heavy particle turbulent dispersion in random flows. 

Equation (|3.6p clearly matches the two limiting behaviours for very small and very large 
separations. In the former case, inertia dominates the small-scale velocity statistics with 
respect to the underlying smooth fluid velocity, and caustics lead to a pure power-law 
behaviour. 



SrV ^ [Snur^S''' ^ V^, 



7(St) 



(3.7) 



where the local Stokes number has attained its dissipative limit St{R) — > St. 

In the latter case, at very large scales R ^ R*{St) inertia is subleading, and the typical 
velocity difference between particles is close to the fluid velocity increment. 



SrVsi 



V^tSRu; 



(3.8) 



At intermediate scales, for large Stokes, St > 1, inertia brings a non-trivial dependency 
via the scale-dependent Stokes number, St{R), and we expect a pseudo power-law scaling: 

SRVst- (Sru)''^^"^^^^ ^ R''^^"^'^^^^/^; ?7<i?<i?*(50. (3.9) 

Summarising, we propose that at changing the Stokes and Reynolds numbers, different 



Turbulent pair dispersion of inertial particles 



9 



R = R (St) 




1 



St 



Figure 4. Sketch of the different regimes expected in the parameter space of inertia, St, and 
scale separation R. The curve St{R*) = 1 separate the region of low inertia St{R) < 1, region 
(C), from the regions where inertia is important St{R) > 1, regions (A) and (B). Further we 
can distinguish the separation regime where inertia is important and particle velocity difference 
larger than the fluid one at the same scale, region (A), from the intermediate regime where 
inertia is still important particle velocity difference is smaller that the corresponding fluid one, 
region (B). Separation between region (A) and (B) is given by the curve SnVst ~ 5ru. For 
relative dispersion of pairs of Stokes St starting at a given separation R, one typically starts 
from the corresponding position in this plane and then evolve upwards along the vertical arrow. 



regimes governing the particle velocity statistics can be distinguished. The relevance 
of such regimes of the particle velocity statistics for the associate relative dispersion 
dynamics can be easily explained with the help of the sketch reported in Fig. [H In the 
parameter space of inertia and scale separation {St, R), we can distinguish three regions 
depending whether inertia is strong or weak, and whether particle velocity difference is 
large or not with respect of the fluid velocity difference at comparable scale. In agreement 
with what commented before, we pose that the curve St{R*) = 1 distinguishes the region 
of weak {St{R) < 1) and strong inertia {St{R) > 1). 

Regime (A) is such that inertia is important since the scale R*{St) ^ rj, and moreover 
the typical particle velocity increments are larger than the fluid increments. In the region 

(B) , inertia is still important but particle velocity increments are depleted with respect to 
the fluid increments. This typically happens for large Stokes numbers, and in our DNS is 
visible only for very large separations R{t) of the highest Stokes St = 70. Finally, regime 

(C) is characterised by a weak inertia, which appears only in the filtering factor for 
the velocity large-scale amplitude and possibly in sub-leading corrections to the tracers 
relative dispersion. 

Even for the largest value of the Reynolds and Stokes numbers achieved in our DNS, it 
is very difficult to disentangle quantitatively the above mentioned regimes, because of the 
closeness of the three relevant scales, 77, R*{St), and L. Still, the quality of the fit shown 
in Fig. [3] using the global functional dependence given by Eqn. (|3.6p makes us confident 
that the main physical features are correctly captured. Before closing this section, we 
note that there is no reason to assume that the functional form entering in the pseudo- 
power law scaling in the inertial range, j{St{R)), in (|3.9p is equal the one characterizing 
the scaling in the viscous range, j{St), in Eqn. (j3.7p . Hint for this observation come 
from results obtained in IBec et al. 20081 for random flows, where a very high statistical 
accuracy can be achieved: there, depending if the underlying fluid velocity is spatially 
smooth or rough, a slightly different functional form has been found. 

The previous analysis gives us a clear quantitative picture of the scale and velocity 
ranges where caustics play a role in the particle dynamics. For example, for moderate 
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Stokes numbers, we have important departure from the tracers statistics only for very 
small scales, i.e. caustics gives a singular contributions to the particle velocity increments 
inside the viscous range; then, at larger scales, the particle velocity scaling become indis- 
tinguishable from the tracer velocities. Clearly, for such Stokes, no important corrections 
for particle separation evolution is expected with respect to the usual Richardson disper- 
sion observed for tracers. This is because particle pairs tends to separate, and very soon 
all pairs will attain separations where their velocities are very close to the underlying 
fluid. On the other hand, for very heavy particles, those with Stokes time falling inside 
the inertial range of fluid velocity statistics, the contribution from the caustics will be felt 
also at relatively large scales, up to i? R*{St). Pair separations attain such scales when 
the initially large relative velocity difference has relaxed and become smaller than the 
corresponding fluid one — crossing from region (A) to (B). Notice that at i? ~ R*{St), 
we have that SR'Vst — V^iSr^u, i.e. there is a non-trivial effect from inertia. Moreover, 
for large Stokes, at scales R > R*{St), particle velocity increments are smaller than the 
fluid counterparts, indicating an important depiction induced by the Stokes drag on the 
particle evolution. 

It is clear from the above discussion that new physics should appear for the value of in- 
ertia and scales separation of region (B). This regime — that we can not access with the 
present data — is the one where a new law of pair separation should appear as recently 
suggested by IFouxon fc Horvai 20081 A discussion of the dispersion regimes of inertial 
particle pairs follows in the next section in terms of the time behaviour of the mean 
square separation distance. 



4. Dispersion regimes and corrections due to inertia 

In this section we analyse the effects of inertia on the mean square separation of heavy 
particle pairs with a given initial separation distance, Rq at time t — to, as a function of 
the Stokes number: 

{{R{t)f I i?o, to)st = - X2{t)f}su (4.1) 

where in the left-hand side the average is performed over all pairs of particles such that 
\Xi(to) — X2(to)\ = Ro- The study of the relative dispersion of small, neutrally buoy- 
ant tracer particles has recently been the subject of renewed interest. This has been 
motivated by the fact that very accurate — highly resolved in time and space — data 
have become available, experimentally (" Ott fc Ma nn 2000*. Bourgoin et al. 2006 ) and nu- 



merically dYeung fc Borgas 2004 , B iferale et al. 200 5, Biferale et al. 2006^ . These studies 
have confirmed what was known since the works of lRichardson 19291 and lBatchelor 1952|, 
i.e. the existence of different dispersive regimes for tracer pairs in turbulent flows, de- 
pending on the value of their initial distance and on the time scale considered. 

When released in a statistically homogeneous and isotropic, turbulent flow with an 
initial separation Rq in the inertial range for fluid velocity, i.e. -q ~ Rq <^ L, tracer pairs 
initially separate according so the so-called Batchclor regime, 

{{R{t)f\Ro,to)st=o-Rl+C{eRof'''t'; t„ « - to) « , (4.2) 

where C is supposed to be a universal constant, and e is the average kinetic energy 
dissipation of the flow. This ballistic regime appears because initially tracers separate as 
i/the underlying velocity field were frozen, and it lasts for a time scale that is a function 
of the initial separation itself, ts = {Rq/e) (see lBatchelor 19521 [Bourgoin et al. 2006D . 
After such a transient initial time, the relative separation dynamics forgets the initial 
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Figure 5. Mean square separation versus time, for heavy particles at changing St and the initial 
distance Rq. Time is normalised with the Kolmogorov time scale r,,. Left panel: St = 0,0.6, 1, 
and 3.3; initial distance iio G [0 : 1]?;, Run I. Error bars due to statistical fluctuations are of the 
order of the symbol size. Notice that the two largest Stokes numbers show a time lag interval 
where separation proceeds faster than tracers. Inset: ratio between the heavy particle separation 
and the tracer data, Q{t) versus time, for St = 0.6, 1, and 3.3. Same symbols as in the body of 
the figure are used. Right panel: mean square separation versus time, but with a larger initial 
distance, Rq £ [4 : 6]rj. Stokes numbers are the same as in the left panel. Notice that now only 
the dispersion of particle pairs with St = 3.3 exhibits a small departure from the underlying 
fluid, as shown by the Q{t) indicator in the inset. For the smaller Stokes, typical size of caustics 
is smaller than the initial separation i?o, and particle pairs therefore separate as fluid tracers 
do. 

conditions and tracers separate explosively with a power law behaviour given by the 
Richardson law: 



mt)Y 



3. 



tB<^it- to) < Tl , 



(4.3) 



^0, to){st=o) ' 

where g is known as the Richardson constant. As set out in Monin & Yaglom fMonin fc Yaglom 2007[ 
the tracer separation PDF — that will be discussed later — has a similar scaling behaviour 
in these ranges. 

A remarkable fact of the Richardson dispersion (|4.3p is the disappearance of the 
dependence on the initial separation i?Oj an effect also dubbed intrinsic stochasticity 



(E & Vanden Eijnden 2000), which is just the signature of the non-Lipschitz nature of 
the velocity field driving the separation between tracers, when their mutual distance is in 
the inertial range of fluid velocity statistics. The experimental and numerical validation 
of the previous prediction (|4.3p has proved to be particularly difficult, the main reason 
being the strong contamination from viscous and large scale effects in the tracers dynam- 
ics. To overcome these problems, a series of techniques have been developed, including 
the study of doubling time statistics; i.e. the probability distribution function of the time 
needed for a pair to double its separation (jBoffetta fc Sokolov 2002) IBiferale et al. 2005|) . 
Thanks to these techniques, a fairly good agreement on the value of the Richardson con- 
stant has been achieved. Here, we want to study how the tracer behaviour is modified by 
the presence of small-scale caustics in particular and by inertia effects in general, for the 
case of heavy particle pairs. Standard direct measurements of the moments of separation 
as a function of time will be considered, while application of doubling time statistics is 
left for future studies. 

In Fig. [51 we show the behaviour for the mean square separation at varying the Stokes 
number, and for two values of the initial separation. We start with data at the lowest 
resolution, i.e. Run I at Re\ ~ 200, and for moderate Stokes numbers, St ^ C(l). Initial 
distances are chosen equal to Rq < rj (left panel) and Rq G [4:6] 77 (right panel). 

If the initial distance is small enough (left panel), the presence of caustics in the particle 
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Figure 6. Mean square separations versus time, for pairs with St = 10, 30, 70, at Rex = 400. 
Left and rigfit panels refer to the two initial distance iio G [0: l]rj and -Rq £ [4:6]77, respectively. 
Error bars due to statistical fluctuations are of the order of the symbols size. Tracers (solid lines) 
are also shown for comparison. Notice the ballistic behaviour for the heavy particle separation 
observed in the caustics dominated time interval. For very large time-lags a Richardson-like 
behaviour starts to develop but with a less intense overall speed of separation, due to the 
depletion effects of the Vg^' prefactor in the particle velocity increments for large Stokes numbers. 
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velocity field at initial time gives a very remarkable departure from the tracer behaviour. 
At increasing the Stokes number, such departure is more and more evident, and it lasts 
for a time lag which becomes longer and longer. For the highest value of the Stokes 
number shown in the left panel of Fig. [5] {St = 3.3), a sensible difference from the tracer 
behaviour is observed over almost two decades: t G [0.1 : 10] t^. A way to better visualise 
the departure from the tracer statistics consists in plotting the mean square separation 
for heavy pairs of different Stokes numbers, normalised to the tracer one, that is 



Qit) = 



St 



(4.4) 



((^W)^)(5*=0)' 

This quantity is shown in the insets of Figure [5l For heavy pairs starting at i?o — ^ and 
with St = 3.3, the relative difference is as large as 10 at its maximum for t Tri- However, 
such effect becomes progressively less important if wc start the separation experiment 
from larger initial distances as shown in the right panel of the same figure. This is because, 
at these same Stokes numbers, the deviation of particle velocity difference with respect 
to the underlying fluid, due to caustics, has already decreased. This is equivalent to state 
that, for these Stokes numbers, the typical size of caustics is smaller than the initial 
separation i?o, and particle pairs therefore separate as fluid tracers do. At larger time 
lags, whatever the value of the initial separation, the Richardson dispersion regime is 
recovered. 

We now consider what happens for larger Stokes numbers. In Fig. [51 we show the results 
for the mean square separation of St = 10, 30, and 70 and for the large Reynolds number. 
Rex ~ 400. Both initial distances. i?o S [0 : 1]?7 and Rq € [4 : Gjry are displayed. As one 
can see, for the large value of St — 70, the tracer-like behaviour is never recovered, 
and even the separation of pairs starting with the largest distance Rq is affected. The 
transient regime dominated by the caustics invades the whole inertia! range: since particle 
pairs need a very long time to decrease their initial velocity difference to the value 
of the fluid increment at the corresponding scale, they separate with a quasi-ballistic 
behaviour: (i?^(t))5t oc t^. 

The above scenario can be interpreted in terms of caustic- dimensions. At any value of the 
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inertia, there exist a spatial length, of the order of the scale R*{St), which identifies the 
typical spatial size of caustics, i.e. the range of scales where particle velocity increments 
are uncorrelated from the underlying fluid velocity field. If the initial pair separation 
Rq is taken inside this region (left panel of Fig. ^ , particle pair separation starts much 
faster than for fluid tracers, because of the much more intense velocity differences felt by 
the pairs inside the caustics. When particle pairs reach a separation larger than R*{t), 
they start to be synchronised with the underlying fluid velocity, recovering the typical 
Richardson dispersion. However, if the initial separation is larger than the caustics size, 
the evolution of inertial particle pairs is almost indistinguishable from the tracers. Finally, 
whether or not a Richardson- like behaviour is recovered for very large inertia, may depend 
on the Reynolds number also. In the limit of larger and larger Reynolds, at fixed Stokes 
number, one may expect a final recovery of the fluid tracers behaviour even for very 
heavy particle pairs. 



4.1. Mean-field approach to heavy particle dispersion 

The turbulent relative dispersion of fluid tracers can be easily modelled by applying K41 
scaling theory to the fluid velocity increments governing particle separation dynamics 
(see, e.g.. lOuellette et al. 2006|) . Indeed, if R{t) is the tracer separation vector at a given 
time, its evolution is completely specified by the equation 

R{t) = u{Xi,t) -u{X2,t) = SRu{R,t), (4.5) 

together with the initial condition i2(io) = -Ro- Hence, we can directly write an equation 
for the root-mean-square separation r{t) = (|i?(t)p | Rq, to)^/"^ 

r = - {R{t) ■ Snu{R{t),t) \ i?o, to) , with r{to) = Rq. (4.6) 
r 

We next assume the following mean-field closure for the right-hand side: 

{R{t) ■ dnu{R{t),t) I i?o, to) « {R^ I Ro, to)'/' {R ■ Sru) = r 5f (r), (4.7) 

where {r) is the first-order Eulerian longitudinal structure function of the underly- 
ing homogeneous and isotropic turbulent flow. According to K41 phenomenology, this 
structure function behaves in the inertial range as (r) ~ Cei/3ri/3, where C IS an 
order-unity constant. This closure finally leads to 

r = Ce^^^r^/^, so that r{t) ^ Ir^^^ + {2C/3) e^/^{t - to)]^^^ . (4.8) 

Such an approximation gives a complete qualitative picture of the time evolution of the 
mean square separation between tracers. In particular, it encompasses the two important 
regimes of relative dispersion: when (t — to) ^ ts = (3/2C) e~^^^RQ^^ , a Taylor expansion 
of the solution (|4.8p gives the Batchelor regime r(t) ~ Rq + C {eRoY^^ (t — to), while 
when (t - to) > ts, one recovers Richardson's law r(t) ~ (2C/3)^/^e^/^ t^/^. 

In the case of inertial particles, the number of degrees of freedom to describe the 
dynamics is obviously increased: the separation between two heavy particles obeys 

R{t) = -— \R{t) - Snu{R,t)] . (4.9) 

In order to derive mean-field equations one has to track simultaneously the average dis- 
tance and velocity difference between particles. For this we follow the same spirit as for 
tracers and introduce the particle velocity structure function v(t) = (|^/jV(t)|^ | i?o, to)"'^''^, 
where 6RV{t) — R.{t) is the velocity difference between the two particles. One can pro- 
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ceed as previously to write from (|4.9p exact equations for r{i) and v{t): 



{R ■ 5ru) 



- - {SrV ■ 6ru) 



(4.10) 
(4.11) 



where for the sake of a lighter notation the indication of conditional ensemble averages 
was dropped. It is worth noticing that the root-mean-square velocity difference v{t) 
evolves with a dynamics that resembles closely that of heavy particles. However, v{t) 
does not coincide with the time derivative of the mean distance r{t). It is thus useful to 
rewrite the above equations introducing a sort of transverse particle velocity component 
w defined as 

f = V -w . (4-12) 
We can write an exact equation also for the evolution of w 



1 , w 1 

w = w — [2v — w) 

r, r T. 



-{R-dRu)--{SRV-6Ru) 

r V 



(4.13) 



Of course, equations (|4.1ip . (|4.12p . and (|4.13p are not closed without supplying the 
correlation between the particle evolution and the underlying fluid. As in the case of 
tracers, the first unclosed term appearing in the right-hand side of (j4.13p is approximated 
by (|4.7p . The next unclosed term involving the correlation between fluid and particle 
velocity differences is approximated by 



{SrV ■ Sru) « {\SrV\Y^^ {\Sru\Y/'' = vSl^^r) 



(4.14) 



where 5*2 (r) denotes the full second-order structure function of the fluid velocity field. 
When r is in the inertial range, K41 phenomenology implies that S'2(r) oc (er)^/'^. Finally 
these approximations lead to a closed set of equations for the time evolution of the average 
separation and velocities r, v, and w 



1 



Ts 

w — (2v 



1/3 



w 

w) — 

r 



'-Be 



1/3^1/3 



(4.15) 
(4.16) 

(4.17) 



where B and C are positive order-unity dimensionless constants, reflecting the lack of 
control on the prefactors of the scaling laws in the closures (|4.7|) and (|4.14|) . This system 
of equations is supplemented by the initial conditions r{to) ~ Ro, v(io) = (I'^flo^l^)^^^ 
and ti;(to) = ^'(^o) ~ {Ro • Sr„Y) / Rq, which clearly depend on the dispersion experiment 
under consideration. It is worth noticing that this system of equations reduces to the 
mean-field equation (|4.8p for tracers in the limit of vanishing inertia —>■ 0. 

Similarly to the case of tracers, the crude approximation (|4.15[) - (|4.17p of the evolution 
of the root-mean-square distance between heavy inertial particles is able to capture the 
main features of the separation time behaviour. In Fig. [71 we show the result of the 
numerical integration of the set of equation (j4.15p - (|4.17p obtained by an appropriate 
choice of the free parameters (see figure caption), together with DNS data from Run II, 
for two different large values of the Stokes number. For fixed initial separation and at 
increasing the intensity of the caustics velocity increments in the initial condition (i.e. at 
increasing inertia), the transient deviation from the Richardson behaviour become more 
and more evident at intermediate times (of the order of the Stokes time Tg, not shown). 
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Figure 7. Mean square separation versus time, for Run II and Stokes numbers St = 10 (left 
panel) and St = 70 (right panel). Data are shown for the available choices of the initial sep- 
arations : f?o < 7?, fio £ [4 : 6] r\, and i?o £ [9 : 11] t\ (from bottom to top). Symbols stand for 
DNS data, while solid lines are mean-field solutions. Adimensional prefactors are _B = and 
C ~ 2. Notice that time is made adimensional with the Stokes time Ts. Vertical dashed lines 
mark t/r^ = 1. 

Clearly such a simple approach can be valid only in a limited region of the phase space, 
where the initial conditions are, at least, at the edge of the inertial range so that K41 
scaling is correct for the fluid velocity second-order increments. Moreover, the matching 
scale where particle velocity increments become of the order of the fluid increments has 
to fall in the inertial range too: if this is not the case, then pairs enter the regime where 
inertia is important but particle relative velocity is small at scales of the dissipative range, 
and the mean-field closure proposed above becomes inadequate. A detailed quantitative 
comparison of the realm of applicability of the mean- field approach — including effects 
of viscous scales and small-scale caustics — will be the object of future work. 



4.2. Cross-over between relaxation to the fluid velocity and Richardson behaviour 

Despite its simplicity, the mean-field approach described above is able to correctly re- 
produce the pair dispersion of heavy particles with initial data in the inertial range. We 
might wonder if one can draw an even simpler qualitative picture of pair dispersion. For 
this, we consider the behaviour of particle pairs with moderately large Stokes numbers, 
for which inertia plays an important role for the initial transient and the Richardson 
behaviour is slowly recovered well inside the inertial range of scales. For simplicity, we 
assume that the scale where the fluid and particle velocity becomes of the same order, 
SrV ~ Sru, and the scale R*{St), where inertia ceases to be important, are very close. 
As it is clear from the sketch of Fig. |4l this may not be always the case because of the 
effect of the normalisation factor Vg^ for large Stokes: in the picture, it corresponds to 
Stokes number with a narrow transient region (B). 

The general picture then goes as follows. Initially particles separate almost ballistically 
during a time which is of the order of (or larger than) the time needed by their initial, 
caustics-dominated, velocities to relax to the fluid velocity. After that time, particles 
behave as tracers and reconcile with a standard Richardson dispersion. This is a first 
order approximation since (?) the fiuid fiow actually correlates to the particle dynamics 
already at very small times, and (ii) inertia effects are present up to large times as 
previously discussed. Nevertheless, such an approximation should give the two correct 
qualitative asymptotic behaviours, at small and large time scales. Since we consider 
moderately large values of the Stokes number, the initial typical particle velocity can 
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Figure 8. Mean square separation versus time from DNS data of Run II, for three different 
values of the Stokes numbers. Notice that time is normalised with the Stokes time Ts. Solid 
lines represent the initial, almost ballistic evolution due to the exponential relaxation of velocity 
statistics, and the dashed lines correspond to the Richardson regime. With a suitable tuning 
of the free parameters, here the Richardson constant g and the initial velocity increment value 
{{SRV{to))'^), both temporal behaviours are reproduced. 

be assumed to be much larger than the fluid velocity, i.e. |(5_rV| ^ Under these 

hypotheses, there is an initial time interval during which difference between particle 
velocities obeys SrV « ~{SrV)/ts [see gUJ], and thus dRV(t) - i6RV{to)) e-(*-*«)/^-. 
As a consequence, the mean square separation between particles evolves initially as: 

ilR'it)] I i?o,to) =Rl + 2Ts{R{to) ■ SRV{to)){l - e-(*-*°)/^0 

+ T^{{6BV{to)f){l - e-(*-*o)/-=)2. (4.18) 

This should be approximately valid up to a time scale, in the inertial range, where 
|<5_rV'| ~ \Sru\ ~ {eRY^^: it is easy to show that such a time scale is proportional to the 
particle response time r^. For larger times, inertia effects become subdominant and heavy 
pair dispersion suddenly gets synchronised to a Richardson like regime. Nevertheless, this 
Richardson regime has started only after the previous relaxation has ended, that is at a 
distance much larger than the original separation Rq of the particle pair. The combination 
of this initial exponential relaxation of heavy particles with moderately large inertia, plus 
the later standard Richardson diffusion are the two main features due to inertia in the 
inertial pair dispersion. This is indeed confirmed by Figure ([S]), where we compare DNS 
data for mean square separation, with the two phenomenological regimes just described, 
for which we have assumed that {R ■ <5_RV(to)} — 0- As we can see, the main qualitative 
trends of the small and large time behaviours are very well captured. 

4.3. Subleading terms in the Riehardson regime 

We have seen in previous subsections that the most noticeable effect of inertia on the 
mean pair dispersion is a long transient regime that takes place before reaching a Richard- 
son explosive separation (|4.3p . and that this regime is due to the relaxation of particle 
velocities to those of the fluid. As we now argue, at larger times — corresponding to 
regime (C) — there is still an effect of particle inertia that can be measured in terms 
of subleading corrections to the Richardson law. To estimate these corrections, let us 
assume that in the mean-flcld equation (|4.16p . the term stemming from the fluid velocity 
(j^i/3^1/3 jg jj-^^nj]^ larger than the inertia term TgV. This is true when St{r) ^ 1, i.e. at 
times t when r{t) ^ R*{St). In this asymptotic, one can infer that the transverse velocity 
component w is much smaller than the total velocity v, so that r ~ u (see eq. (|4.12|) ). In 
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Figure 9. Large-time behaviour of the mean square separation normalised to that of tracers as 
defined from Eqn. (|4.4|) for Run I and various values of the Stokes number as labeled. Deviations 
from the fluid tracer Richardson law behave as {t/Ts)~^ ■ The limiting value Q{T) has been chosen 
different from unity as effects of inertia are still present at the largest scale of the flow. 

the spirit of the weak inertia expansion derived in |Maxey 1987| we next write a Taylor 
expansion of (|4.16p to obtain 

f « « « Ce^'\^'^ - r,(|(d/dt) 5,w|2)i/2 « Ce^/^ri/^ - T,(|<S,a|2)i/2^ (4.I9) 

where 5j.a — dr{dtU + u ■ Vit) denotes the increment of the fluid acceleration over the 
separation r. Next we assume scaling invariancc of the turbulent acceleration field, that 
is, according to dimensional arguments of K41 theory, \Sra\ ^ ^2/3^,-1/3 Equation (|4.19p 
can then be rewritten as 

r = C£i/3ri/3 (1 - ^ r, £i/3^-2/3^ ^ Ce^/^r^^^ (1 - A St{r)) , (4.20) 

where A is an order-unity constant. The initial condition is given by r{to) = tq where 
the initial separation has to be chosen such such that St{ro) ^ 1. We can next integrate 
the approximate dynamics perturbatively in terms of the small parameter St{ro) by 
expanding the separation as r{t) = po{t) + Pi{t) + p2{t) + ■ ■ ■ ■ The leading order is po{t) = 
[tq^^ + (2C/3) e^^^t]'^^'^ and corresponds to the relative dispersion of a pair of tracers. The 
first-order correction reads pi(t) = —TsE^I'^A ln(po(i)/?'o) Pq^^Y times much larger 
than the Batchelor time associated to the initial separation rg, i.e. for t ^ e~^/'^rg^'^, the 
leading term follows the Richardson explosive law p^it) ~ (2C/3)^/^e^/^t'^/^. This finally 
implies that in the asymptotics t 3> e'^^'^r^^'^ ^ r,, one can write 

r^it)(xgt^ [1- D{t/rsy' In [t/Ts)], (4.21) 

where g is the Richardson constant introduced in Sj4l ^-^d D is an order-unity factor, 
which a priori does not depend neither on the particle Stokes number, nor on the initial 
particles separation. 

This behaviour is confirmed numerically as can be seen from Fig. [9] that gives the 
behaviours at large times of the ratio Q{t) between the mean square separation of heavy 
particles and that of tracers as defined by ()4.4|) . One can clearly see that data almost 
collapse on a line oc 1/t confirming the behaviour (|4.2ip predicted above. Only results 
from Run I arc displayed here. The reason is that the very large time statistics of tracer 
dispersion in Run II is not as well statistically converged, leading to more noisy data. 
The qualitative picture is however very similar. 

To conclude this section, let us stress that we have assumed above K41 scaling to hold 
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for the acceleration field (and thus for the pressure gradient). However it is well known 
that the scaling properties of pressure are still unclear: they might depend on the turbu- 
lent flow Reynolds number and/or on the type of flow (see, e.g., |Gotoh fc Fukayama 2001[ 
IXu et al. 2007p . As stated in lBec et al. 20071 rather than being dominated by K41 scal- 
ing, numerically estimated pressure increments of Run I {Re\ ~ 200) seem to be ruled by 
sweeping, so that \5ra\ u^ms e^l'^r~'^l'^ . One can easily check that this difference in scal- 
ing leads to a behaviour similar to (|4.2ip . except that this time logarithmic corrections 
are absent, and that the non-dimensional constant D depends on the Reynolds number 
of the flow. The present numerical data do not allow to distinguish between these two 
possible behaviors. 



5. Probability density function of inertial particle separation 

We now discuss the shape of the probability density function for both light and heavy 
inertial particles. We focus on the time and scale behaviour of the non-stationary PDF 

VstA^AR^M). (5.1) 

defined as the probability to find a pair of inertial particles (S't, /3), with separation i? at 
time t, given their initial separation Bq at time t^. The case of tracers (St = 0, /3 = 1) 
has been widely studied in the past, either experimentally, numerically and theoreti- 
cally for two and three dimensional turbulent flows (see lRichardson 1929|,IBatchelor 1952|, 

lJuUien et al. 1999[lBofetta fc Sokolov 2002|[Biferale et al. 2005l|Bourgoin et al. 20061ISalazar fc Collins 2009p . 
Following the celebrated ideas of Richardson, phenomenological modelling in terms of a 
diffusion equation for the PDF of pair separation leads to the well-known non-Gaussian 
distribution. 



T'st=o.i3=i{R,t) oc ^ exp 



(5.2) 



(£^3^)- 

which is valid for times within the inertial range r,, <^ t <^ T^, and is obtained assuming 
a small enough initial separation and statistical homogeneity and isotropy of the three- 
dimensional turbulent flow. Here, A is a normalization constant. This prediction is based 
on the simple assumption that, for inertial range distances, tracers undergo a diffusion 
dynamics with an effective, self-similar, turbulent diffusivity K{B) cx BSru ~ £^/^ B'^^^. 
Moreover, it relies on the phenomenological assumption that tracers separate in a short- 
time correlated velocity field. Indeed, it is only if the latter is true, that the diffusion 
equation for the pair separation becomes exact (see IFalkovich et al. 200ip . 

As mentioned before such a scenario may be strongly contaminated by particle inertia. 
The main modifications are expected to be due to the presence of small-scale caustics 
for small-to-large Stokes numbers, and to preferential concentration. Caustics make the 
small scale velocity field not differentiable and not self-similar, as if inertial particles 
were separating in a rough velocity field whose exponent were depending on distance. 
Preferential concentration, leading to inhomogeneous spatial distribution of particles, 
manifests itself as a sort of effective compressibility in the particle velocity field. 

There exists a series of stochastic toy models for Lagrangian motion of particles in 
incompressible/compressible velocity fields, where the statistics of pair separation can 
be addressed analytically. Among these, the so-called Kraichnan ensemble models, where 
tracer particles move in a compressible, short-time correlated, homogeneous and isotropic 
velocity field, with Gaussian spatial correlations (we refer the reader to the review 
IFalkovich et al. 20011 for a description of this model) . It is useful for the sequel to re- 
call two main results obtained for relative dispersion in a Kraichnan compressible flow. 
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We denote with p the velocity field compressibility dcgrcc[3, and with < ^ < 2 the scal- 
ing exponent of the two-point velocity correlation function at the scale r, in rf-dimensions: 
{[u,ir) ~ u,(0)]K-(r) - uj{0)]) ^ Gir«[(rf - 1 + ^ - pO% + ^ipd~ l)nr,/r% For parti- 
clcs moving in such flows, it is possible to show that the pair separation PDF for tracer 
particles follows a Richardson-like behaviour: 



V^,^{R,t) K 



R 



D2 



exp 



-A 



i?2 



t 



(5.3) 



Here pi = p£,{d + / {1 + pCj and D2 — d — fi is the correlation dimension, characterizing 
the fractal spatial distribution of particles. 

A different distribution emerges when the d— dimensional Kraichan flow is differentiable, 
i.e. for ^ = 2; in such log-normal PDF is expected: 



Pf_,,^{R,t\RQ,tQ) oc ^ exp 
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(logiR/Ro) - Xjt - to))- 
2A(t-<o) 



(5.4) 



with A = 2Gi{d -!)(! + 2p) and A = G'i(d - l)(d - 4p). It is worth noticing that in 
the latter case, since the flow is differentiable, the large-time PDF depends on the initial 
data. 

The problem of inertial particle separation in a real turbulent flow presents some 
similarities with the previous toy cases but also important differences. 
First, the effective degree of compressibility — due to preferential concentration of inertial 
particles — , is properly defined only in the dissipative range of scales. For r ^ 77, it is 



equal to the correlation dimension D2 defined as p{r) 



where p{r) is the probability 



to find two particles at distance smaller than r, with r <^ 77. As it has been numerically 
shown in IBec et al. 2007| ICalzavarini et al. 20081 for three-dimensional turbulent flows, 
the correlation dimension depends only on the degree of inertia [St, /3), while it does not 
seem to depend on the Reynolds number of the flow. For r ^ rj, the effective degree of 
compressibility is no longer constant, but varies with the scale. 

Second, the underlying velocity field exhibits spatial and temporal correlations that are 
much more complex than in a Gaussian short-correlated field. Such correlations lead to 
non-trivial overlaps between particle dynamics and the carrying flow topology. As a result, 
it is not possible to simply translate the analytical findings obtained in the compressible 
Kraichnan ensemble to the case of inertial particles: we may expect however that in some 
limits the compressible Kraichnan results should give the leading behaviour also for the 
case studied here of inertial particles in real turbulent flows. 

With this purpose, we first notice that the separation probability density function 
that is valid in the rough case (|5.3p has an asymptotic stretched-exponential decay that 
is independent on the compressibility degree. This suggests that inertial particle PDF 
(|5.ip must recover the Richardson behaviour (|5.2p of tracers in the limit of large scales 
and large times. Coherently with what discussed in previous sections, for large times and 
for scales larger than Rg^, we expect that the heavy pairs (in the limit /? ~ 0) PDF 
recovers a tracer like distribution : 



VstfiiR.t) ~ exp 



-A 



R^R 



Sf 



(5.5) 



For pairs of light particles, there is no straightforward formulation of such a prediction: 



f The compressibility degree p is defined as the ratio p = /S'^ , where oc ((V • u)^) and 
o oc {{S/uf), and varies between p = for incompressible flows, and p = 1 for potential flows. 
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Figure 10. Separation probability density function, P(7?, t|7?o, to), for heavy pairs with different 
Stokes numbers, at changing time. Initial distance is taken Ro £ [3 : 4] 77 for St = (top-left), 
St = 0.6 (top-right), and St = 3.3 (bottom-left) of Run I, and equal to iJo £ [4:6] r; for St = 70 
(bottom-right) of Run II. The related initial distributions are pictorically depicted with a grey 
area. Times shown are: {t — to)/T,, = 1, 6, 18, 36 for Run I and (t — to)/T,, = 1, 6, 18, 36, 86 for 
Run II. 



as we shall sec in the sequel, preferential concentration effects have a strong fingerprint 
on the separation PDF even at large times and large scales. 

In the opposite limit of very small separations, i.e. R rj, one can correctly assume 
that the effective degree of compressibility is constant and therefore apply either the 
small-scale limit for rough flows (|5.3p . or that for smooth flows (|5.4p . depending on the 
scaling properties of the particle velocity field entailed in the value of the exponent 'y{St), 
defined from (|3.4p and related to the caustics. We thus expect 

VstAR. t\Ro. to) ~ R^~'G{t), if jiSt) ^ 1, 

VstAR^ t\Ro, to) - R'^^-^Fit), if j{St) = 1. (5.6) 

Here F and G are two different decaying functions of time t, whose expression can be 
easily derived from (|5.3p - (|5.4p . Notice that for the smooth case, i.e. the small-scale limit of 
the log- normal distribution (|5.4p . we get for the spatial dependency a factor £'2/2 instead 
of the factor D2 of the rough case. This will matter in the case of light particle separation, 
where, due to strong preferential concentration, the probability of finding pairs at a very 
small distances is large enough to allow for a detailed test of the prediction (|5.6p . The 
case of light particles will be discussed in !j6l while we now turn to a discussion of the 
above scenario in the case of heavy particle pairs. 
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Figure 11. Comparison of PDFs at fixed times with data of Fig. (|10|l . Left: early stage of the 
separation process, t — to = Trj- Inertia does not affect small Stokes, St = 0.6 while its effect 
is detectable for St = Z and St = 70. Right: PDFs comparison at a later time, t — to = 36r,,. 
Now the PDF shows some deviations from the tracer behaviour only for St = 70. On the right 
panel the solid line gives the Richardson shape (|5.2|> . Initial separation and Reynolds numbers 
are the same as for Fig. (|10p . The inset shows the PDF evolution for St = 70 at three times, 
(t - to)/T-„ = 36,82, 130. 



5.1. Probability density function of heavy particle relative separation 

We start by analyzing the qualitative evolution of Vst{R, t) at changing time, for different 
Stokes numbers and in the limit /3 = 0. The four panels of Fig. [TO] show the evolution 
of the PDF at different times for pairs with St ^ (tracers), and for heavy particles 
with St = 1, 3.3, and 70. Initially, at t = tQ, all selected pairs are separated by the same 
distance {Rq £ [3 : 6] 77) ; this initial distribution is represented in each figure by a grey area. 
As time elapses, particle separate and reach different scales, depending on their inertia. 
Qualitatively, the PDF evolution is very similar for all moderate Stokes numbers, and the 
PDFs at different moderate Stokes numbers become more and more similar with time. 
However in the case of St = 70 — for which the associate Stokes time Ts falls well inside 
the inertial range — , the PDF shows a long exponential tail for intermediate separation, 
which tends to persist at all observed times. To better appreciate such differences, in Fig. 
[TT] we show the comparison between the different PDFs corresponding to various Stokes 
numbers for two different times: at the beginning of the separation process, {t — to) = r,,, 
and at a later time, (t — to) = 36t^. As one can see, it is only at early times that 
the PDFs for moderate-to-large Stokes, St = 3, 70 differ in a sensible way from the 
tracers. In particular, one can clearly see that many pairs have separations much larger 
and much smaller than those observed for tracers or for heavy pairs with small Stokes 
numbers. The right tails, describing pairs that are very far apart, are just the signature 
of the scrambling effect of caustics. Such strong events are not captured by second order 
moments of the separation statistics that we discussed before, while they clearly affect 
higher-order moments. The left tails, associated to pairs much closer than tracers, are 
possibly due to particles that separate at a slower rate than tracers because of preferential 
concentration induced by inertia. 

Later in the evolution, for (t — tg) — 36t^, only the separation PDF for St — 70 still 
shows important departure from the tracer case; for all the other Stokes numbers shown, 
pairs have had enough time to forget their initial distribution and have practically relaxed 
on the typical Richardson-like distribution. In the inset, we also show the persistence in 
the exponential behaviour for the PDF at St = 70, by superposing the shapes measured 
at three times during the particles separation. 

With the present data the small scale asymptotic behaviour (|5.6p cannot be validated 
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Figure 12. Time evolution of the probability density function of heavy particle relative lon- 
gitudinal velocity, Wst{v,t), during the separation process. Data refers to four different cases: 
tracers pairs St = 0, and heavy pairs St = 1,3, 70, starting with initial distance f?o £ [4 : 6)77. 
PDFs are measured at times {t — to) = [0, 1, 38, 61]t7,, for Run II. Notice the presence of intense 
velocity fluctuations for moderate-to-strong inertia, St = 3, 70, observable at the early stage of 
the separation process. These are the legacy of the caustics distribution. 

for heavy particles. This is due to the limited statistics: very soon after the initial time 
to, there are almost no pairs left with separations R rj. 



5.2. Probability density funetion of heavy partiele relative velocities 

At moderate to large Stokes numbers the separation process of heavy particle pairs is 
largely influenced by the presence of large velocity differences at small scales, that is by 
the presence of caustics in the particles velocity field. In iJ31 we have studied station- 
ary statistics (only first-order moment) of velocity differences between heavy particles at 
changing the distance between particles and their inertia. However it is also informative to 
look at the non-stationary, time-dependent distribution of velocity differences, and more 
particularly to its distribution measured along heavy pairs separation. The relative veloc- 
ity 6iiV{t) — Xi{t) — X2{t) can be decomposed into the projection along the separation 
vector, and two transveral components, here equivalent since the system is statistically 
isotropic. For tracer particles, the statistics of relative velocity and the alignment proper- 
ties of (5fl V'(t) and R{t) have been discussed extensively (see e.g. |Yeung fc Borgas 2004p . 
Here, we focus on the PDF of the relative longitudinal velocity only, which we denote by 
Wst{v,t), where v{t) = [Xi{t) - X2{t)] ■ R{t). For pairs of tracers {St = 0), the initial 
longitudinal velocity distribution is nothing else than the PDF of Eulerian longitudinal 
velocity increments measured at the distance Rq. For pairs of inertial particles, this ini- 
tial PDF clearly coincides with the stationary distribution of velocity differences between 
particles that are at a distance R = |Xi(to) — X2{to)\ G [Rq-.Ro + dRo]. Such a distribu- 
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Figure 13. Time evolution of the mean square sej^aration for two different families of light par- 
ticles {St = 1.2, /3 = 3) and {St = 4.1, /3 = 3). The case of tracers is also shown for comparison. 
Notice that the strong small-scale clustering does not affect the long-time behaviour, except 
through a very small asymptotic slow down. Inset: ratio between the mean square separation 
for light pairs and that of tracers. 

tion has the signature of two mechanisms: (i) at small Stokes numbers, only preferential 
concentration matters and particles probe only a sub-set of all possible fluid velocity 
fluctuations; {ii) at large Stokes numbers, particles arc homogeneously distributed but 
with a velocity field which may be strongly different from the underlying fluid velocity. 

For what concerns heavy pairs, the first effect has not an important signature on small- 
scale quantities. However the second effect clearly becomes visible for moderate to high 
inertia as shown in Fig. 1121 Here we report the longitudinal velocity distributions for 
pairs with initial distance i?o G [4 : 6] r/, and with St = 0, 1, 3.3, and 70; the Reynolds 
number of the underlying flow is Re\ ^ 400. Each panel contains the PDFs measured at 
different times spanning all turbulent timescalcs. At t = io, the importance of caustics is 
manifest for the two largest Stokes numbers, leading to fat tails towards both small and 
large velocity differences. Interestingly enough, the left tail of Wst{v,t), which describes 
approaching events of particle relative motion, is immediately dumped already at {t — 
to) ^ T,^: at the same time, however, the right tail continues to be quite fat for the two 
largest Stokes numbers under consideration. At later stages of the separation process, 
the tendency of large-Stokes-number pairs to wash out approaching events becomes even 
stronger. Indeed, at time {t — to) = 38 t^, the small velocity increments tail has almost 
disappeared for pairs with St = 70. It is worth noticing that at those times (i.e. also at 
those typical scales) , heavy particle velocity differences have already started to be smaller 
than the tracer velocity increments: the larger is the Stokes number, the less pronounced 
are the PDF tails. 

Summarising, because of the different effects of inertia, we observe a very complex 
evolution for the longitudinal relative velocity fluctuations along the trajectories of heavy 
particle pairs. This is certainly a key issue to be considered for stochastic modelling: 
here, as in a standard kinetic problem, both particle positions and velocities need to be 
modelled to quantitatively control the relative dispersion process. 

6. Relative dispersion for light particles 

So far we have considered the relative motion of very heavy particle pairs, for which 
the density contrast (3 with the underlying fluid is zero. In this section wc present results 
on light particles dynamics as described by (|2.2p . for different possible choices of the 
parameters {St,P). 

We discuss how the strong effect of preferential concentration — typically observed in 
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the case of light particles in turbulent flows — might influence the intermediate and long- 
time behaviour of pair separations. In three-dimensional turbulent flows, as we consider 
here, light particles associated to different values of {St, (3) have been observed to always 
possess a positive largest Lyapunov exponent ICalzavarini et al. 20081 this implies that 
light pairs always separate in 3d real turbulent flows. We recall, however, that this is not 
always true: for instance, in smooth two-dimensional random flows, there are values of 
{St, (3) for which the largest Lyapunov exponent can become negative and particles form 
pointwise clusters fsee lBec 2003[) . 

Light particles with moderate inertia and high density ratio (order-unity St, and 
/3 = 3), initially tend to separate much slower than heavy particles with the similar 
Stokes number: this is evident from the much smaller values of the Lyapunov exponents 
measured for light particles, with respect to those measured for heavy particles with 
equivalent Stokes numbers but /3 = 0. Moreover, finite-time Lyapunov exponents show 
large fluctuations, indicating that there are pairs that do not separate even at long times. 
Results on this issue will be reported elsewhere. Clearly, pairs of light particles that do 
not separate do not influence the mean square distances: hence, we do not expect, and 
indeed do not measure, any large differences for the long-time behaviour of {\R{t)\^) st,t3 
for light particles, with respect to the heavy case (see Fig. [T3l) . 

It is natural to ask if light particle strong preferential concentration affects high-order 
moments of relative separation of two initially close particles, and particularly the left tail 
of the separation PDF. Figure [14] shows the time evolution of the separation probability 
density function, 'Pst.i3{R,t\Rn,ti)). Data refer to a case with very intense preferential 
concentration effects and minor influence of caustics {St = 1.2,(3 = 3), and a case 
with milder inhomogeneities in the spatial distribution {St = 0.3, /3 = 2). The initial 
separation PDF was chosen in both cases by selecting particle pairs with initial distance 
Rq G [4:6] 77. A remarkable observation is a strong tendency to flU small separations. In 
other words there are many pairs that reduce their mutual distance even for a very long 
times. The development of the left tail for the strong clustering case {St = 1.2, /3 = 3) is 
consistent with the estimate given by the long-time, small-scale asymptotic expansion of 
the log-normal distribution (|5.6p . V{R, t) ^ R^'^l^~^ , as shown by the straight line in the 
plot. This is the confirmation that the small-scale dynamics of the highly clustered light 
particles evolves as that of tracers moving in a smooth, compressible flow (characterised 
by the same D-i). 

We also remark that if there is high spatial preferential concentration, caustics cannot 
be important. This may have important consequences for the estimation of collision 
kernel of light particles. The approaching events, shown by the left tail in Figure [141 ^.re 
clearly due to the preferential concentration inside vortex-like structures, typical of light 
particles. The right panel show a different case, where preferential concentration is less 
important, leading to a correlation dimension Z?2 = 2. Of course, also in this latter case, 
there are events with approaching pairs, but these become less and less probable with 
time. 

The importance of preferential concentration can also be appreciated by looking at the 
PDFs of longitudinal velocity differences between light particles during the separation. 
We show such distributions for one of the pair family considered above, and we compare 
them with those of the tracers (sec Fig. [T5)) . The important difference between the two 
cases stems from the highly peaked nature of the relative velocity PDF for the strong 
clustered light particle case. The presence of many pairs with almost vanishing velocity 
differences is the signature of a coherent bunch of pairs moving on a strongly clustered 
set. 
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Figure 14. Time evolution of the separation probability density functions of light particles. 
Left panel; St = 0.3, [3 = 2, corresponding to a case where preferential concentration is not 
very effective {D2 = 2). As time elapses, one observes a self-similar filling toward smaller sep- 
arations, in agreement with (|5.6|l (dashed line). Right panel: St = 1.2, /3 = 3, corresponding 
to light particles with strong clustering properties (correlation dimension D2 = 0.8). Again the 
self-similar filling of small scales is consistent with the prediction D2/2 — 1 as depicted by the 
dashed straight line. 
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Figure 15. Probability density functions of relative longitudinal velocity, yVst,fs(v,t), for light 
pairs with St = 1.2,13 = 3, and for tracers {St = 0,/3 = 1). PDFs are measured along the 
separation process at two different times: left panel refers to the initial time, t = to; while right 
panel refers to t — to = 38r^ . 

7. Conclusions 

We have studied the relative dispersion of inertial particles in homogeneous and isotropic 
turbulence from two DNS at resolutions 512'^ and 2048^, corresponding to Rex ^ 200 
and Re\ ^ 400, respectively. We have analysed both heavy and light particle statistics at 
changing the Stokes numbers. We have studied the evolution of mean separations and the 
whole PDFs' shape, both for particle distance and velocity increments at changing time 
and for different typical initial distances. The main results that we have discussed can be 
summarised as follows. Separations of very heavy particles, with Stokes times falling in 
the inertial range of the underlying fluid, are strongly affected by the presence of caustics 
up to times, when the distance between particles reaches scales that are large enough 
for the separation dynamics to be again dominated by the underlying flow velocity. As 
a consequence, strong transient departure from the Richardson diffusion, with a faster 
ballistic regime, is observed. A statistical closure of the equation of motions for heavy 
particle separation is also developed. This model is able to reproduce the main numerical 
findings. 

For light particles, at high density ratio, we observe strong small-scale clustering prop- 
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erties, leading to a considerable fraction of pairs that do not separate at all — although 
the maximum Lyapunov exponent remains positive. In such a case, the non-stationary 
spatial concentration at small scales tends to be higher than the analogous case but with 
a stationary distribution of particles. Such numerical findings open the way to exper- 
imental verifications and gives input to the community involved in modelling inertial 
particle diffusion in applied configurations. 
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